Transcriptional changes in specific subsets of Drosophila neurons following inhibition of the serotonin transporter

The transcriptional effects of SSRIs and other serotonergic drugs remain unclear, in part due to the heterogeneity of postsynaptic cells, which may respond differently to changes in serotonergic signaling. Relatively simple model systems such as Drosophila afford more tractable microcircuits in which to investigate these changes in specific cell types. Here, we focus on the mushroom body, an insect brain structure heavily innervated by serotonin and comprised of multiple different but related subtypes of Kenyon cells. We use fluorescence activated cell sorting of Kenyon cells, followed by either or bulk or single cell RNA sequencing to explore the transcriptomic response of these cells to SERT inhibition. We compared the effects of two different Drosophila Serotonin Transporter (dSERT) mutant alleles as well as feeding the SSRI citalapram to adult flies. We find that the genetic architecture associated with one of the mutants contributed to significant artefactual changes in expression. Comparison of differential expression caused by loss of SERT during development versus aged, adult flies, suggests that changes in serotonergic signaling may have relatively stronger effects during development, consistent with behavioral studies in mice. Overall, our experiments revealed limited transcriptomic changes in Kenyon cells, but suggest that different subtypes may respond differently to SERT loss-of-function. Further work exploring the effects of SERT loss-of-function in other Drosophila circuits may be used help to elucidate how SSRIs differentially affect a variety of different neuronal subtypes both during development and in adults.

serotonergic signaling. Relatively simple model systems such as Drosophila afford more 23 tractable microcircuits in which to investigate these changes in specific cell types. Here, we 24 focus on the mushroom body, an insect brain structure heavily innervated by serotonin and 25 comprised of multiple different but related subtypes of Kenyon cells. We use fluorescence 26 activated cell sorting of Kenyon cells, followed by either or bulk or single cell RNA sequencing to 27 explore the transcriptomic response of these cells to SERT inhibition. We compared the effects 28 of two different Drosophila Serotonin Transporter (dSERT) mutant alleles as well as feeding the 29 SSRI citalapram to adult flies. We find that the genetic architecture associated with one of the 30 mutants contributed to significant artefactual changes in expression. Comparison of differential 31 expression caused by loss of SERT during development versus aged, adult flies, suggests that 32 particular serotonergic neurons. This, coupled with the genetic tools available in flies, affords a 66 technically tractable platform for the molecular interrogation of serotonergic circuits and in 67 particular, specific subsets of post-synaptic neurons that receive serotonergic inputs. proposed to regulate different behavioral outputs [35][36][37]. For experiments involving drug-induced SERT blockade (Fig. 5), female flies were sorted on the 88 day of eclosion and maintained on 1% agar + 5% sucrose + 1% blue food dye, with or without 89 The final set included 93084 SNPs, which were transformed into heterozygous variants for the 161 demultiplexing of F1 samples (i.e. alleles were modified from 1/1 to 1/0). The same VCF file was 162 used for demultiplexing of all experiments. The genotypes that were not used in a particular 163 experiment/sample were used as negative controls. Raw sequencing reads and the VCF file for 164 demultiplexing will be available at the NCBI repository (upload to GEO in progress). 165 Single-cell data analysis was performed using Seurat (v4.1.1) [49,50]. Single-cell transcriptomes 166 were filtered using the following criteria: (1) transcript count ³ 1000; (2) maximum percentage of 167 mitochondrial transcripts £ 20%; (3) we also removed cells that were classified by demuxlet as 168 "doublets/ambiguous", and cells that were assigned to the genotypes that were not used in the 169 given experiment.  GPCRs (Dh44-R1, Proc-R, CCHa2-R, Ir76a) (Fig. 1C, D and Supp. Table T1). When genes 214 were plotted by chromosomal position, however, there was a striking concentration of DEGs on 215 the same arm of the 2 nd chromosome (chr2R) as the dSERT 16 DNA lesion (Fig 1E). Drosophila 216 have only 3 chromosomes that house most of their genome, and some of these observations 217 may represent true findings. However, the buildup on chr2R suggests that at least some of the 218 observations may derive from disruption of genomic DNA rather than changes in serotonergic 219 signaling. 220 221 Though SMART-seq libraries feature increased sensitivity to lowly-expressed transcripts, they 222 necessitate pooling of RNA from all cell-types within the collected population and may result in 223 washout of cell-type specific changes. To investigate the transcriptomics of each KC subtype 224 independently, we followed a recent single cell RNA-seq strategy in which all samples and 225 replicates are pooled and processed together [38,44]. We generated dSERT 16 and dSERT 4 fly 226 lines with GFP expressed in KCs as above, but included an additional element unique to each 227 biological replicate: a 3 rd chromosome derived from independent WT strains available from the 228 Drosophila Genetics Research Panel (DGRP) [45]. Because transcripts derived from DGRP 229 chromosomes bear SNPs, single cells can be bio-informatically traced to genotype-of-origin 230 post-hoc ( Fig. 2A). This allowed us to pool all replicates of both control and mutant samples for 231 dissociation, FACS, library prep, and sequencing, thereby minimizing long-standing issues of 232 technical variability between individual replicates that contribute to bias in RNA-seq data. 233 Dimensionality reduction (Supp. Fig. S1) resulted in robust clusters for two sub-populations for 234 KCα/β (KC_AB1, KC_AB2), two for KCγ (KC_G1, KC_G2), and one for KCα'/β' (KC_ABp1) (Fig.  235 2B). Running pseudobulk differential expression between mutant and control cells collapsed by 236 cell-type revealed 33 significant changes. Some changes were cell-type specific (e.g. SK in 237 KC_G1 and CG31690 in KC_AB1), and many were observed in multiple cell-types (e.g. prom, 238 Cbp53E, CG42392, Pgant9) (Fig. 2C, D and Supp. Table T2). For those DEGs that were 239 identified as cell-type specific such as SK, we detected robust transcript expression in most of 240 the clusters, lending credence to the hypothesis that the DE observed is in fact specific to a 241 particular cell-type (Supp. Fig. S2). When visualized in pseudo-Manhattan plots (Fig. 2E), 242 however, the bias of DEGs to chr2R was even more pronounced than for SMART-seq (Fig. 1E), 243 highlighting their possible artefactual provenance. The DEGs on chr2R appear to lie in two 244 positional "columns" -one ~7.5 Mb away from dSERT, and one that is immediately adjacent to 245 the dSERT 16 deletion. One of the DEGs immediately adjacent to the deletion is an eye-specific 246 gene, prom, that is not expressed in WT KCs. By extension, we concluded that upregulation of 247 the prom transcript in dSERT 16 is likely to represent an artefact caused by the deletion of 248 regulatory DNA adjacent to dSERT and prom. 249

250
To explore the possibility that more precise mutations in dSERT might be less disruptive and 251 generate fewer artefactual hits, we generated a new mutant allele using CRISPR [67] to 252 precisely excise ~2.6kb DNA coding for most of the first and second transmembrane domains 253 and simultaneously induce a frameshift in the CDS. We reasoned that even if the resultant 254 mRNA could code for a partial dSERT protein, it would be topologically inverted in the plasma 255 membrane (Fig. 3A). Fly lines bearing the deletion, termed dSERT TMKO , were outcrossed six 256 times to w 1118 . The presence of the deletion was confirmed by PCR-sanger sequencing, and 257 behaviorally in that this line phenocopies the sleep deficit found in dSERT 16 (data not shown). 258 We then built fly lines as in the previous experiment, using the new dSERT TMKO allele and 259 second chromosomes derived from w 1118 as controls, in place of dSERT 16 and dSERT 4 , 260 respectively. Sample prep, scRNA-seq, and data processing (  Table T3). However, in this dataset there 263 is no pronounced enrichment of DEGs on chr2R (Fig. 3F). Importantly, some of the DEGs on 264 chr2R in the previous (dSERT 16 ) dataset, including those immediately adjacent to dSERT, such 265 as prom, are absent from this dSERT TMKO dataset (Supp.  Table T4). Interestingly, some genes (e.g. LysRS) were shared with the previous (day 0) 279 dataset, while Cbp53E, a gene identified in the dSERT 16 day 0 dataset but not found in the 280 dSERT TMKO day 0, reappeared in this dSERT TMKO day 4-6 dataset.  Table T5). As predicted, there was no "pileup" of these observations on 291 chr2R (Fig. 5E). We hypothesized that loss of dSERT activity during both development and adulthood, rather 382 than development alone, might further alter the DE profile. To test this, we repeated the scRNA-383 seq protocol using flies that had been aged for 4-6 days rather than freshly-eclosed (day 0). We 384 again observed some changes across multiple cell types (i.e. LysRS, CG42260), as well as 385 some that were cell-type specific. Among these, the cell surface recognition molecules beat-IIa 386 and side DE in KC_G2 could, similarly to dpr1 in KC_AB1 in the experiment with day 0 flies, 387 represent homeostatic changes to maintain proper connectivity. However, the total number of 388 DEGs seen in the aged flies was similar to that seen with newly eclosed flies. 389

390
To further explore the effects of dSERT inhibition in the adult, we fed WT flies the SSRI 391 citalopram (CIT) or vehicle (VEH) for 4-6 days and repeated our scRNA-seq workflow. We 392 uncovered a new set of DEGs, most of which were observed only in the major KCα/β subtype 393 (KC_AB1) and which did not show significant overlap with those detected using mutants. It is 394 possible that off-target effects of CIT dominate these observations, and drug specificity may be 395 tested in future experiments by feeding CIT to dSERT mutants. It is also possible that the 396 decrease in SERT activity caused by citalopram was less pronounced than the complete block

Conflicts of interests 502
There are no competing financial interests in relation to the work described. showing log2 fold-change (L2FC) for DEGs in dSERT 16 versus dSERT 4 at day 0, comparing bulk 844 sequencing (Fig. 1) and the initial scRNA-seq data (Fig. 2) analyzed using "pseudobulk" to 845 collapse all clusters into one artificial "cell-type" for comparison with the bulk dataset.  Table T1. DE table for  end that includes a non-coding exon and upstream regulatory DNA. The dSERT4 genetic backgroundmatched control contains a 278 bp deletion but does not signi cantly alter protein expression or behavior compared to WT [62]. B) Sample preparation for bulk sequencing. Flies contained the Mef2(P247)-gal4 driver and UAS-nls.GFP marker for expression in KCs, and were homozygous for either dSERT16 (mutant) or dSERT4 (control) on the second chromosome. Flies were dissected and pooled by genotype, then dissociated and FACS-sorted in parallel to select for GFP-labeled KCs, followed by isolation of RNA for bulk RNA-seq (SMART-seq). C) Volcano plot showing differential expression between dSERT16 and dSERT4 groups. DE genes include those encoding the transcription factors Lim1 and Achi, the channels Ork1 and Ppk29, the GPCRs Dh44-R1, Proc-R, CCHa2-R, and Ir76a, the calcium binding protein Cbp53E, and genes implicated in neuronal development (Trim9, Mis12)  scRNA-seq of KCs from dSERT16 and dSERT4 ies, in immediately-eclosed (day 0) ies.
A) Flies used for scRNA-seq contained one of six unique 3rd chromosomes derived from different DGRP wild-type lines, as well as the markers and dSERT alleles used for bulk seq. Two and four different DGRP lines per group (dSERT16 or dSERT4 , respectively) were created and served as biological replicates.
Brains from all lines were dissected, pooled, and dissociated together, then FACS-sorted to select KCs used for scRNA-seq. B) t-SNE dimensional reduction showing distribution of cells in this dataset among transcriptionally-de ned clusters (see methods) representing KCg cells (KC_G1, KC_G2), KCa/b (KC_AB1, KC_AB2), and KCa'/b' (KC_ABp1). C) Volcano plot from "pseudobulk" analysis (by cluster) of DEGs between dSERT16 and dSERT4 . Observations are color-coded (as in B) by the KC-type in which they were identi ed. D) Venn Diagram showing overlap of DEGs identi ed in the major cell clusters. Cbp53E, CG42392, and CG33143 were identi ed as DE in multiple cell types. E) DEGs plotted by chromosomal locus as in Figure 1E. A skewed localization of DEGs to chr2R is notable. Figure 3 dSERTTMKO scRNA-seq, in immediately-eclosed (day 0) ies.
A) Cartoon depicts the independently-derived dSERTTMKO deletion compared to dSERT16 . B) Flies used for this scRNA-seq experiment were homozygous for dSERTTMKO or a WT dSERT allele derived from w1118 and expressed with the same transgenes for isolation of KC cells as in Fig 1 and 2. Each y was marked by a different DGRP 3rd chromosome variant, and t-SNE plot shows the color-coded distribution of cells by KC cell-type as in Fig. 2. C) Volcano plot as in Fig. 2C from "pseudobulk" analysis (by cluster) of DEGs between mutant and control. Observations are color Drosophila serotonin coded (as in B) by the KC-type in which they were identi ed. D) Venn Diagram showing overlap of DEGs identi ed in the major cell clusters. CG42392 and LysRS were identi ed as DE in both KC_G1 and KC_G2. E) DEGs plotted by chromosomal position as in Fig. 2F. In contrast to Fig. 2, observations are not concentrated on chr2R. Figure 4 scRNA-seq for dSERTTMKO vs controls in aged (day 4-6) ies.
A) Flies harboring the dSERTTMKO or WT dSERT alleles (in the control line w1118 ) were aged for 4-6 days then processed for scRNA-seq as in Figure 3  scRNA-seq in aged ies treated with an SSRI.
A) Flies with WT dSERT alleles were treated with citalopram (CIT) to block SERT protein activity or vehicle alone (VEH). Each y contained one copy of 2nd and 3rd chromosomes derived from a unique DGRP line and transgenes for marking KCs as in previous gs. B) t-SNE plot indicating the distribution of cells by cell-type. C) Volcano plot from "pseudobulk" analysis (by cluster) of DEGs between mutant and control.
Cell-type speci c DEGs include Lgr1 and Ddc in KC_AB1 and Hsp26 and Hsp70Bc in KC_G2, none of which were identi ed in previous experiments. D) Venn Diagram showing that there is no overlap of DEGs identi ed in the major cell clusters. E) DEGs plotted by chromosomal position as in previous gs. Similar to Figs. 3 and 4, observations are not concentrated on chr2R.

Figure 6
Correlation of genes identi ed as DE between datasets.
A) Correlation plot showing log2 fold-change (L2FC) for DEGs in dSERT16 versus dSERT4 at day 0, comparing bulk sequencing (Fig. 1) and the initial scRNA-seq data (Fig. 2) analyzed using "pseudobulk" to collapse all clusters into one arti cial "cell-type" for comparison with the bulk dataset. Concordant genes signi cant in both datasets are plotted in a larger font, and colored purple. Genes signi cant in only the bulk or scRNA-seq datasets are colored red or blue, respectively. Diagonal dark grey dashed line represents 1:1 correlation between datasets. The lighter grey horizontal and vertical lines represent 1.5 fold-change cutoffs for genes of interest. B) Correlation plot between dSERT16 and dSERTTMKO day 0 scRNA-seq datasets. Genes are color-coded by KC type as in previous gures. Genes not signi cant in either dataset are plotted with reduced opacity. Genes signi cant in at least one dataset are plotted with normal opacity. While there are many genes with L2FC of the same sign in both datasets, most are only signi cant in one dataset (smaller labeled points). C) Correlation plot comparing data derived from newly eclosed (day 0) vs aged ies (day 4-6) using the dSERTTMKO . Genes are plotted as in B. CG42392 and LysRS in KC_G1 were signi cant in both datasets (larger labels and points), with DE in the same direction (downregulated). D) Correlation plot between aged dSERTTMKO (d4-6) and aged ies fed citalopram (CIT). One gene (Hsp26) was DE in both datasets, although in a different cell type in each dataset and therefore not highlighted.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.